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Abstract: The near-side ridge observed in A+A collisions at RHIC has been described as 
arising from the radial flow of Glasma flux tubes formed at very early times in the collisions. 
We investigate the viability of this scenario by performing a non-perturbative numerical 
computation of double inclusive gluon production in the Glasma. Our results support the 
conjecture that the range of transverse color screening of correlations determining the size 
of the flux tubes is a semi-hard scale, albeit with non-trivial structure. We discuss our 
results in the context of ridge correlations in the RHIC heavy ion experiments. 
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1. Introduction 

At very high energies, the collision of two nuclei can be conveniently described as the 
collision of two Color Glass Condensates (CGCs) [jy, [2], [J. In this description, the wave- 
functions of the nuclei are comprised of high occupation number classical gluon fields at 
small x coupled to static light cone color sources at large x. The separation between fields 
and sources evolves with energy; one obtains evolution equations for multi-parton corre- 
lations in the nuclear wavefunctions called JIMWLK renormalization group equations (2[. 
After the collision, the coherence of the nuclear wavefunction is lost; large x sources are no 
longer static and become sources for multi-particle production in the forward light cone § . 
For early times r ~ l/Q s , where Q s is the saturation scale, the non-equilibrium produced 
matter with high occupation number of order l/a s has remarkable properties g, |6[ and is 
termed the Glasma [Qj. For reviews, see Ref. ||. 
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Recently, high energy factorization theorems were derived for inclusive multi-gluon 
production in a rapidity interval Ay < l/a B 
expressed very compactly as p 



1C] in A+A collisions. The result can be 



d n N n 



d 3 P! • • • d 3 p n 



dN 



[D Pl ][D P2 ] Z y [ Pl ]Z y [p 2 ] -g- 



dN 



d 3 P„ 



(1.1) 



The Z's are gauge invariant weight functionals that describe the distribution of color 
sources at the rapidity of interest. They are obtained in full generality by evolving the 
JIMWLK equations from an initial rapidity close to the beam rapidity. In the large A c limit, 
the weight functionals Z can instead be obtained from the simpler mean field Balitsky- 
Kovchegov (BK) |llj equation. In this case, they can be represented as non-local Gaussian 



distributions in the sources [12]. For large nuclei, without significant small x evolution, one 
recovers the local Gaussian distributions of the McLerran-Venugopalan (MV) model Q. 
We emphasize that the validity of eq. ([O]) is restricted to the kinematics where all the 
produced particles are measured within a rapidity interval < l/a s from each other, so that 
we can evaluate Z at this same rapidity y\ « • • • ~ y n « y. The generalization to larger 
rapidity separations is non-trivial because one needs to account for the gluon radiation 
in the region between the tagged gluons. A formalism describing arbitrarily long range 



rapidity separations has been developed recently [13]. For simplicity, we shall not consider 
it further here. 

The leading order single particle distributions in eq. ( |1.1|) are given by 



dN 



d 3 p 





PUP2 


LO 





lim / d 3 x d 3 y (Sg - iE p )(d° y + iE p ) 

x J24(PVx(P) A% cL [pi,P2](x) Ar l [Pi,P2](y) ■ (1.2) 



For each configuration of sources p\ and p 2 respectively, of each of the nuclei, one can 
solve classical Yang-Mills equations to compute the gauge fields A c ^[pi, p 2 ] in the forward 



light cone 0, H When our expression for the corresponding single inclusive 

distribution is substituted in eq. (|1 . 1| ) , and the distributions are averaged over with the 
distributions Z, one has determined from first principles (to all orders in perturbation 
theory and to leading logarithmic accuracy 1 in x), the n-gluon inclusive distribution in 
high energy A+A collisions at proper times t ~ l/Q s - As noted previously, eq. (1.1) is 
valid only for AY < l/a s ~ 3 — 5 units of rapidity in A+A collisions at RHIC and the 
LHC respectively. 

In the nuclei, before the collision, the typical range of color correlations is of the order 
of the inverse saturation scale Q$ l , where Q s _1 < Aqcd -1 - The saturation scale at a given 
transverse position in the nucleus depends on the two dimensional transverse projection 



of the nuclear matter distribution. In eq. (Fl), it appears in the initial conditions for the 
Z functionals; the energy evolution of the saturation scale is determined by the JIMWLK 



Note that what we mean by leading log of x here means resumming the leading dependence in rapidity 
between the observed gluons and the projectiles; not between the different produced gluons. 
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renormalization group equations. Because the saturation scales in the two nuclei are the 
only scales in the problem (besides nuclear radii), their properties determine the energy 
and centrality dependence of the inclusive observables. 

The expression in eq. (|0|) is remarkable because it suggests that in a single event- 
corresponding to a particular configuration of color sources-the leading contribution is from 
n tagged gluons that are produced independently. The coherence in n-gluon emission is 
generated by averaging over color sources that vary from event to event. Because the range 
of color correlations in the transverse plane is of order 1/Q S , the high energy factorization 
formalism suggests an intuitive picture of correlated multiparticle production arising from 
event by event fluctuations in particle production from ~ R 2 A /1/Ql = R\Qg color flux 
tubes of size 1/Q S . 

The color flux tubes, called Glasma flux tubes, are approximately boost invariant in 
rapidity because of the underlying boost invariance of the classical fields 2 . Besides providing 
the underlying geometrical structure for long range rapidity correlations, the Glasma flux 
tubes also carry topological charge |J, which may result in observable metastable CP- 
violating domains |l9|]. While eqs. Ql.l| ) and (1.2) describe particle production from the 



Glasma flux tubes, they do not describe the subsequent final state interactions of the 
produced gluons which become important for times r 3> 1/Qg- 

If, as widely believed, the produced matter thermalizes by final state interactions, 
these will not significantly alter long range rapidity correlations because the effects of 
fragmentation, hadronization or resonance decays are typically restricted to Ay ~ 1 -C 
l/a s . However, radial flow will have a significant effect on the observed angular correlations. 
This is because particles produced isotropically in a given flux tube will be correlated by 
the radial outward hydrodynamic flow of the flux tubes. Ideas on the angular collimation 



of particle distributions by flow were discussed previously in the literature [20, 21 1. When 



combined with the long range rapidity correlations provided by the Glasma flux tubes, they 



provide a natural explanation [22, 23] of "ridge" like structures [^4|, £5|] seen in the nearside 



spectrum of two particle correlated hadron pairs in A+A collisions. A similar structure is 
seen in three hadron correlations as well [26]. 

These structures were first seen in near side events with prominent jet like structures, 
where the spectrum of associated particles is observed to be collimated in the azimuthal 
separation A<p relative to the jet and shows a nearly constant amplitude in the strength 



of the pseudo-rapidity correlation Arj up to Ar] ~ 1.5 [24]. The name "ridge" follows 
from the visual appearance of these structures as an extended mountain ridge in the Arj- 
Aip plane associated with a narrow jet peak. This collimated correlation persists up to 
Arj ~ 4 ]P7t] . An important feature of ridge correlations is that the above described 
structure is seen in two particle correlations without a jet trigger and persists without 
significant modification for the triggered events [pf| . These events include all hadrons 
with momenta px > 150 MeV. In such events, a sharp rise in the amplitude of the ridge 



is seen [28, 29 1 in going from peripheral to central collisions. A number of theoretical 



For rapidity separations > l/a s , one expects significant violations of boost invariance from quantum 
corrections [tL3l. 
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models have been put forth to explain these ridge correlations [21, 30]. For a recent critical 
evaluation of these models, we refer the reader to Ref. [|3l| ]. 

In this paper, our primary purpose is to discuss how the Glasma flux tube picture 



arises ab initio from solving eqs. ( |1 . 1[ ) and ( |1.2| ). In Ref. [22], it was argued that the two 
particle correlation function 

c _ / d 2 N 2 \ / dN \ / diV \ 

2 \<iy p d 2 p T dy q d 2 q T / \ dy p d 2 p T J \ dy q d 2 q T J ' 

in the Glasma flux tube picture took on the simple geometrical form 

C2(p,g) 1 n ^ 

/ di v W div \ K2 5 ± Q2' U - 4j 

\ dy p d z ip T I \ dy q d J q T / 

where the right hand side of this relation is simply a non-perturbative constant k 2 mul- 
tiplied by the ratio of the transverse area of the flux tube to the transverse overlap area 
of the nuclei. This nice geometrical identity however is a conjecture based on a perturba- 
tive computation whose regime of validity is the kinematic region pt,Qt 3> Qs, where one 
expects additional contributions, besides those arising from strong sources, to contribute 
significantly. Similar arguments, based on perturbative computations were used to gener- 
alize the result in eq. ( |1.4[ ) to 3-particle correlations [32] and subsequently, even n-particle 



correlations [33|. In the latter paper, it was shown that the general structure of n-gluon 
emissions is a negative binomial distribution. 

Because the geometrical structure of the correlations is so striking and with manifest 
consequences, it is important to establish that is has a validity beyond perturbation theory. 
The perturbative expressions extended to lower momenta are strongly infrared divergent 
so it is not obvious whether they are regulated in that regime by the confinement scale of 
the order of Aqcd or instead a semi-hard scale of order Q^ 1 arising from the weak coupling 
albeit non-perturbative dynamics of the Glasma fields. We will tackle the problem head-on 
in this paper and investigate the form of eq. (|1.3| ) in this paper. We will do so by solving 
the Yang-Mills equations in eq. (|1.2| ) in full generality numerically on a space-time lattice. 
Our non-perturbative results from the lattice simulations are valid 3 in the entire kinematic 
domain of transverse momenta-from the smallest infrared scale given by the lattice size to 
the largest ultraviolet scale given by the lattice spacing. In physical units, these correspond 
to the inverse nuclear size and large momenta (pt,Qt S> Q s )- 

One can of course, on dimensional grounds, always express the r.h.s of eq. ( |1.4| ) as 
shown. However, k 2 in general can depend on the dimensionless ratios Q s /pt> Qs/qt, the 
relative azimuthal angle Aip between the two gluons, and the dimensionless combinations 
Q s Ra and m/Q s , where m is an infrared scale of order Aqcd and Ra is the nuclear radius. 
We will examine the dependence of n 2 on these parameters carefully and discuss what they 
tell us about the structure of correlations. We find that the Glasma flux tube picture is 
valid; however, while the size of transverse correlations is still a semi-hard scale, it is not 



3 We emphasize however that at very large momenta there will be additional contributions to the distri- 
butions beyond the purely classical one discussed here. 
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simply Q s and does display some sensitivity to infrared physics. Our results are important 
for quantitative explanations of the ridge as resulting from Glasma flux tubes. In the semi- 
quantitative computation of Ref. [23], K2 was implicitly taken to be a free parameter and 
fit to the data-we will compare our result to this value. We will also compare our result 
to the value obtained by comparing the Glasma multiplicity distribution [33] to PHENIX 
data on multiplicity distributions at RHIC p4|]. 

In section ^, we will review perturbative computations of multi-particle production in 
the Glasma. The non-perturbative computation will be discussed in Section ||. We will 
briefly outline the numerical approach, clearly stating the lattice parameters, their relation 
to physical parameters, and the approach to the continuum limit in the transverse and 
longitudinal co-ordinates. We will then systematically study the dependence of our results 
on the saturation scale, the nuclear size and the infrared cutoff m. In section ^, we will 
discuss the physical implications of our results and further refinements for quantitative 
comparison with experiment. We will end with a brief summary of our key results. Some 
details of the numerical computation are discusssed in an appendix. 



2. Perturbative results for multi-particle production 



The perturbative computation of two, three and ra-gluon correlations in the Glasma have 
been discussed elsewhere [22, E32L 53]. For completeness, we will briefly review the results 



here. The correlated two particle inclusive distribution can be expressed as 



C(p,q) 



1 



4(2vr) e 



£ (\MTkp,q)\ 2 )-(\M a x ( P )\ 2 )(\M a ;(q)\ 



(2.1) 



a,a';A,A' 



where the classical contribution to the amplitude for the production of a pair of gluons 
with momenta p and q is 



M&,(p,q) = e>K {q) P 2 q 2 A^{p)A»' a (q) , 
M a x (p) = e x Jp)p 2 A^ a (p). 



(2.2) 



Here the e's are the polarization vectors of the gluons and a, a' are the color indices of the 
gauge fields. The average (• • • ) in eq. (2.1) is an average over the color configurations of 
the two nuclei; this average will be discussed further shortly. 

The gauge fields have a very non-trivial, non-linear dependence on the sources pi and 
p2, which evolves as a function of the proper time r. As mentioned previously, they can be 
determined by numerically solving Yang-Mills equations for r > with initial conditions 
at r = given by the gauge fields of each of the nuclei before the collision [14, 17]. 
We will discuss these numerical solutions further in the next section. For large transverse 
momenta pr 3> Q s , however, the equations of motion can be linearized and one can express 
the classical gauge fields produced in the nuclear collision as H |35|, 36, [37j 



p 2 A^ a (p) = -if abc g 3 



d 2 k 7 



L"(p,k T , 



p\(^t)p c 2 {vt - k r ) 



•4 (PT 



(2.3) 
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Here f a f, c are the SU(3) structure constants, is the well known 4 Lipatov vertex and pi,p 2 
are respectively the Fourier transforms of the color charge densities in the two nuclei ||. 
From eq. (1.1), the average in eq. ( |2.lD corresponds to 



(O) = J [Dp 1 Dp 2 ]Z[p 1 ]Z[p 2 ]0[p 1 ,p 2 ] . (2.4) 



In the MV model 



Z[p] = exp f- J d 



V*^| , (2.5) 



where p can be either p\ or p 2 . As we will discuss further in the next section, the color charge 
squared per unit area p? can be expressed simply in terms of Q s . We will consider this 
Gaussian model in the rest of this paper 5 . For these Gaussian correlations, in momentum 
space, 

[p a (k T )p b (kir)) = (2tt) 2 p 2 a 5 ab 8(k T - k^) . (2.6) 
Using eq. (^3|), fl2j) and (gj) in eq. @), one obtains H, 

16 tt< g^Qi p^q^ 



It is instructive to express the result in eq. (2.7) in terms of the inclusive single gluon 



spectrum. This result, due originally to Gunion and Bertsch [BO], has been been recovered 



previously in the CGC framework M, 35, 41, 03 and is known to have the form 



dN \_S ± ( ff Vj 4 N C (N? - 1) ^ /pf 



dy p d 2 PT / 4vr 4 g 2 p\ \Q, 

Substituting eq. on the r.h.s of eq. (^7]), one obtains 



which is the result we noted in eq. (|1.4|), with 6 



Identifying the theoretical error on k 2 is difficult at this stage from the perturbative com- 
putation. 

Computing k 2 non-perturbatively by numerically solving the Yang-Mills equations is 
the primary objective of this paper. This quantity is not a pure number but can contain 



4 The components of this four vector are given explicitly by L + (p, kr) = — jpr, L (p, kr) = 

5 In the simplest treatment of small x evolution, based on the Balitsky-Kovchegov equation |nj, Z[p] 
can also be modelled by a Gaussian, albeit one with a non-local variance fj9[ . 

6 In the computation of Ref. [H, there were numerical errors which gave a significantly larger perturbative 
estimate for «2. 
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interesting structure. In the perturbative calculation, K2 is independent of the relative angle 
A<p at very large transverse momenta pr, q? 3> Q s . However, even in the perturbative 
calculations, one notices that there are finite azimuthal correlations between gluons as one 
goes away from the limit of asymptotically large transverse momenta. It is important to 
understand these correlations at momenta of interest to experiment to ascertain whether 
they have any phenomenological significance. Further, K2 can depend non-trivially on the 
transverse momenta of the pairs and on the energy and centrality of the nuclear collision. 
Finally, we want to determine whether K2 is infrared finite. This was the case for the 



analogous factor for the non-perturbative single gluon distribution at any finite time [14]; 
there is no guarantee that this should be the case for multi-gluon distributions. 

Before we end this section, we should mention that perturbative computations were 



also performed for 3-gluon [32] and n-gluon correlations [33]. The result can be nicely 



summarized in terms of the cumulants of the multiplicity distribution as [33] 



d n N 



dyx d 2 pTi ■■■ dy n d 2 p Tn 



(n- 1)! 



fcn-l 



dN 



dyi d 2 p Tl 



x • • • x 



dN 



dy n d 2 p T „, 



where 



k = C 



(iV 2 - 1)Q 2 S S ± 
2vr 



(2.11) 



(2.12) 



Here £, is a non-perturbative coefficient 7 to be determined by numerical solutions of Yang- 



Mills equations. The expression in eq. (2.11) corresponds to a negative binomial distribu- 



tion. One can extract £ from fits to the multiplicity distribution data [34] and compare our 
non-perturbative result for K2 to this value thereby testing the validity of eq. ( [2,llD for the 
2nd cumulant. Our approach can be straightforwardly be extended to higher cumulants 
but the computations are numerically intensive and will not be considered further here. 



3. Non-perturbative computation 

The two gluon cumulant, on dimensional grounds, can always be expressed in the form 
given in eq. fl2.9|) . However, as we have noted previously, the dimensionless coefficient K2 
can in general be a function of Q s /pt, Qs/qt, A</j, Q s Ra and m/Q s , where m is an infrared 
regulator scale that can be added to the MV model. If we assume that m ~ Aqcd, the high 
energy asymptotics of our formalism is m/Q s <C 1 and Q s Ra 00. For A+A collisions 
at realistic RHIC and LHC energies, one expects m/Q s ~ 0.2 — 0.1 and Q s Ra = 35 — 70 
respectively. We will return to the discussion of the infrared scale m later at the end of 
sec. |. 

For asymptotically large px, qr S> Q s , we anticipate, on the basis of our perturbative 
results, that K2 — > constant. For p?, qx < Q s , we will determine the dependence of K2 on 
these ratios. In particular, we would like to determine whether the double inclusive gluon 
distribution is rendered infrared safe by the non-linearities of the Yang-Mills fields that 
may generate a mass scale (analogous to a plasmon mass) for finite times r > 1/Q S - 

7 In the perturbative computation for the two and three particle correlations, this number was taken to 
be C = 2. 
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3.1 Numerical approach 

The numerical solutions of the classical Yang-Mills equations have been discussed at length 
elsewhere 14, 15, |l7]] and we will only outline the approach followed here for completeness. 



In the MV model the Yang-Mills equations 

[D fM ,F^] = r. (3.1) 

have the source currents 

J^ = ^ + p 1 (x r ,x-)+^- /02 (x r ,x+), (3.2) 

where the color charge densities p\$, of the two nuclei are independent static color sources 
localized on the light cone pi i 2(xr,x =F ) ~ pi^xt)^^). They are distributed with the 
Gaussian probability distribution 

(p a (x T )p\y T )) = g 2 p 2 A 5 ab 5 2 (x T - y T ), (3.3) 

where and yy are vectors in the transverse plane. The initial conditions for the solutions 
of these equations in the forward light cone in the Fock-Schwinger gauge (A T = 0) are given 
by 

A i \ T=0 = A\ + Al (3.4) 
A\ =0 = l f[A\,Ai), 

with 

4,2 = -^1,2*9^2 and V^A lj2 = -gp h2 . (3.5) 

where 

Vi, 2 = V T e iAl ' 2 . (3.6) 

It is convenient in the forward light cone to use the co-ordinate system (r, rj, xy), where 
r 2 = t 2 — z 2 and r\ = 0.51n((i + z)/(t — z)). The Yang-Mills equations in this co-ordinate 
system are manifestly boost invariant with the gauge fields independent of r], namely, 
A^ = #(r, x^). The symbol Vzp denotes path ordering in the =F directions of the Wilson 
line V\ t 2 respectively, which for nucleus 1 (nucleus 2) have an implicit dependence on x~ 
(x + ) in the solution. We will return to this point at the end of the subsection. 

A subtle point we would like to emphasize is that the Wilson lines V, in the gauge field 
solution given by eq. (^5j) for nucleus 1 (nucleus 2) before the collision, are path ordered in 
x~ (x + ). This feature of the continuum solution is implemented by constructing the Wilson 
lines in consecutive steps representing a discretization of the longitudinal coordinate. On 
each lattice site x^ one constructs random color charges with a local Gaussian distribution 

(fiWftW) = 5 a W(x T - y T )^ , (3.7) 

with the indices k, I = l,...,N y representing a discretized longitudinal coordinate. Nu- 
merical calculations of the single inclusive gluon distributions previously used = 1, 



-8- 



whereas the derivation of the analytical expression in eq. ( |3,5|) required an extended source 
corresponding to the limit N y — > oo. 



Our normalization is chosen to give 
52 (pg(x r )rf (y T ) 



5 ab 5 2 (x T -y T )g'^ A . 



2 ..2 



(3.8) 



The Wilson lines for the initial condition are then constructed from the sources eq. ( |3.7| ) 
by solving a Poisson equation and exponentiating it to give 



Vi,2(xt) 




pI' 2 ( x t) 
V| + m 2 



(3.9) 



Here we have introduced an additional infrared regulator m into the MV model; some kind 
of infrared cutoff is required for inverting the Laplace operator. The single inclusive gluon 
distribution, integrated over transverse momenta, turns out to be weakly dependent on the 
infrared cutoff. One of the central aims of this paper is to study whether this is also the 
case for correlations. The parameter m plays an important role in the interpretation of our 
results and we will return to it in our discussion of the results of our computations. We will 
return to the physical interpretation of the role of the parameter m later in our discussion. 



For large N y the charge densities pk in eq. (|3.7| ) become small, and the individual elements 



in the product (3.9) approach identity. This is precisely the procedure that leads in the 



N„ 



oo limit to the continuum path ordered exponential in eq. (p7 



The previous numerical computations of the single inclusive spectrum defined the 
gluon multiplicity in a manner that slightly deviates from the definition based on the 
LSZ formalism in eq. (|1.2|). In |i"7fl , taking advantage of the equipartition of energy in the 
classical theory, only the electric field parts of the numerical solution were used. Assuming 
a free dispersion relation the gluon multiplicity was taken to be 



dN 



1 



rTr 



^(kr^-kr) + tE^{\l t )E^(-\lt) 

T 



(3.10) 



d 2 k T dy (2tt) 2 |k T | 

where the electric fields E l = tA{ and E v = A^jr are time derivatives of the gauge 
potentials, defined here with the explicit factors of r coming from the formulation of the 
theory in r, 77-coordinates. In Refs. [|i~4| , [D| the numerically obtained E- and ^4-fields were 
used to determine the dispersion relation of the interacting theory. This was then in turn 
used to determine the multiplicity, unlike the free dispersion relation assumed in eq. (|1.2|). 
In this work we use a definition that corresponds to eq. (p^) and take the gluon spectrum 
as 

1 



dN 



1 



d 2 k r dy (2tt) 



rTr 



rk 



-E l {k T )E l (-k T ) + r|k T |^(k T )^(-k T ) 



+ -^r^(k T )^(-k T ) + i^A,(kT)A,(-k T ) 



+ i 



+ i 



E i {\z r )A i (--k T ) - Ai(k T )E\-k T ) 



E*(k T )A v (-k T ) 



A v (k T )E^(- 



(3.11) 

(3.12) 
(3.13) 
(3.14) 
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where the fields are evaluated in the 2-dimensional Coulomb gauge. The interference 
terms ( |3.13| ) and (3.14) are odd under the transformation kj- — > — ky. They do not 



contribute when the gluon spectrum is averaged over the angle of kp, integrated over 
k^ or averaged over configurations of the sources. Thus neglecting them was justified 
when one was interested in the single inclusive gluon spectrum. In the case of two-gluon 
correlations they cannot, however, be neglected 8 . Due to these interference terms, the 
symmetry n(kr) = n(— kr) does not hold configuration by configuration, but only on the 
average. Thus the correlation function C(pt,^t) ^ C(pt,— qr)- In particular, there is a 
peak in the correlation at p-r = qr, which without the interference terms ( 3.13 ) and ( 3.14| ) 



would, by symmetry, imply a similar peak at px = — qr- The main numerical effect of 
including these terms is that the backward peak at p-r = — qr is significantly depleted. 

3.2 Parameters in the computation 

From the discussion in the previous subsection, the parameters in the numerical lattice 
computation are 

• g 2 /J,A, the root mean square color charge density. 

• N y , the number of points in the discretization of the longitudinal (x~ or rapidity) 
direction. 

• The infrared regulator m. When m = 0, the Poisson equation is solved by leaving 
out the zero transverse momentum mode. This procedure corresponds to an infrared 
cutoff given by the size of the system. 

• The lattice spacing a. 

• The number of transverse lattice sites Nt, giving the transverse size of the lattice 
L = N T a. 

Of the dimensionful parameters a, g 2 \iA an d m, only the dimensionless combinations g 2 HAO- 
and ma appear in the numerical calculation. The continuum limit a — > is taken by 
letting Nt — > oo such that g 2 HAa — > and g 2 fiAL = g 2 ^LAdNT is a constant. This constant 
is determined by the physics input of the calculation. We have t^R\ = L 2 and, as we 
shall discuss in the next subsection, g 2 [iA is simply related to the saturation scale Q s . 
The physical combination Q s Ra relevant at RHIC and LHC energies will translate into 
corresponding values of g 2 fiL in our computation. In sec. |J we will explicitly translate 
our numerical results into physically relevant numbers. We note that there is another 
dimensionless combination m/g ,2 ^ J 4-we will study the sensitivity of our results to this ratio. 

3.3 Relating g 2 [iA and Q s on the lattice 

The root mean squared color charge density g 2 [iA that appears in our Glasma computations 
can be simply related to the saturation scale Q s extracted from deeply inelastic scattering 
experiments. In these experiments, the cross section for a virtual photon scattering off 



3 We thank F. Gelis for pointing this out to us. 
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a high energy hadron or nucleus is simply related to the dipole cross section of a quark- 
antiquark pair scattering off the hadron. This "dipole" cross-section is determined, in the 
CGC formalism, by the correlator of two Wilson lines in the fundamental representation 
of each of the nuclei before the collision as 



C F (x r - y T ) = (Tr W(xrMyr)> 



(3.15) 



with the expectation value () evaluated with the distribution of the sources. For Gaussian 
correlators in the MV model of sources extended in the longitudinal direction, 



( P a ( XT ,x-)p b (y T ,y~))=gW(x T -y T )5(x- -y-)^ A (x-), 

the Wilson line correlators can be computed analytically up to a logarithmic infrared cutoff 
that must be introduced in solving the Poisson equation in eq. (|3.5[). The result is 



C F (x T ) m N c exp ( — xxfdn(m|x T |) 



with 



X = 9 



dx fi\(x ). 



(3.16) 
(3.17) 



Alternately, the saturation scale Q s , in the fundamental representation, is defined as the 
momentum corresponding to the maximum of k^C* (ky), where C F (kx) is the Fourier 
transform of eq. ( 3.15| ). In this manner, one can relate the saturation scale to the root 
mean square color charge density. 

The saturation scale defined previously 
is an inverse correlation length associated 
with the correlator of two Wilson lines in 
the fundamental representation. In a nuclear 
collision, it is the saturation scale in the ad- 
joint representation that is relevant. It is 
defined, equivalently as the momentum cor- 
responding to the maximum of k^C^kr), 
where C A (kx) is the Fourier transform of 
the correlator of two adjoint Wilson lines 




C A (x T ) = {V ab (x T )V ab (Q)) 



(3.18) 



Figure 1: Adjoint representation Wilson line 



with 



correlator for different 
k T /Q s . From Ref. ||. 



N v as a function of 



K 6 (x T ) = 2Tr \t a V\-x T )t b V{* T ) 



(3.19) 



With some algebra, in the large N c limit, the adjoint representation correlator can be 
related to a higher correlator of fundamental representation Wilson lines 



C A (x T ) 



Tr 



y f (x T )y(o) 



i 



(3.20) 
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In the Gaussian MV model, one obtains 

C A (x T ) ps (N 2 - 1) exp ( — xx^ln(m|x T | 



Nr. 



8tt' 



(3.21) 



The saturation scale Q s in the adjoint representation, in identical fashion to the fundamen- 
tal case, is defined as the momentum corresponding to the maximum of ky(> (kr), where 
C (kr) is the Fourier transform of eq. ( 3.2C| ). The adjoint Wilson line correlator is plotted 
in fig. [I] The two saturation scales are approximately related by the ratios of the Casimirs 
of the representations, namely, Q 2 
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Figure 2: Left: Q s /g Ha versus N y . From Ref. 43 . Right: Dependence of the ratio Q s /g g, on 
m/Q s for a fixed Q s — 1 GeV and N y = 50. These are the conversion factors used to obtain the 
m-dependence of our results for fixed Q a in fig. || 

Our definition of the saturation scale is not sensitive to the exact shape of the correlator 
for very large or small transverse momenta, and for a Gaussian correlator it reproduces the 

~ 2 

GBW saturation scale as 1/Rq = Q s ■ For small momenta, the fundamental and adjoint 
correlators look like Gaussians, which is the form used in the "GBW" fit of DIS data in 
Refs. [Q. For large momenta there is a power law tail 1/ky that differs from the original 
GBW fits, but resembles more closely the form required to match smoothly to DGLAP 
evolution for large Q 2 pEBT. 



In Ref. HU, the relation between Q s and g 2 ^A was studied extensively employing 
numerical lattice techniques. Discretizing the longitudinal distribution of sources as de- 
scribed in eq. (3/7) and constructing the Wilson line as in eq. ( |3.9| ), it was shown that 
Q s ~ 0. 57 g 2 fi a for N y = 1. This value was used in the numerical Glasma simulations of 
the single inclusive gluon spectrum. As we shall see in the following the gluon spectrum, 
when expressed in terms of the scaling variable Pt/Qs, is to a very good approximation 
independent of N y . This is however not the case for the double inclusive spectrum, which 
will have a stronger dependence on N y . 

Due to the m-dependence of the single inclusive spectrum, the ratio n>2 defined by 
eq. ( |1.4[ ) will depend on m. To study this m-dependence for larger values of N y (closer to 
the continuum limit in the longitudinal direction), we will need to invert the computation 
of Q s as a function of g 2 [x to obtain, as a function of m/Q s , the required values of g 2 \i 
corresponding to a fixed value of Q s . The result of this exercise for Q s = 1 GeV and 3 GeV 
and using the numerical method employed in [43] is shown in fig. ||. 



- 12 - 



4. Results 



We begin the discussion of our numerical results with the single particle spectrum because 
we are interested in the ratio of the double inclusive spectrum to the single particle spectrum 
squared. Some of the results for the single particle spectrum are new and have not been 
published elsewhere. We will then proceed to discuss the double inclusive spectrum and 
state our results for K2 defined in eq. Ql.4|) , 



4.1 Single inclusive spectrum 

The single inclusive multiplicity in A+A collisions was computed extensively previously by 
solving Yang-Mills equations numerically. We note that the single inclusive multiplicity 
computed gives excellent fits to the RHIC multiplicity data |l6|, with values of Q s for 
gold nuclei that agree to ~ 10% with those extracted from extrapolations of HERA data 



to RHIC [ 46 1 . In this sub-section, we will address some details of the computation that 
have not been presented previously and are relevant for the computation of k-i- 

In Ref. KJ, the dependence of the single inclusive gluon distribution as a function of 



N y was not computed. We have done it here and the result is shown in fig. y. It shows that 
the single inclusive spectrum, scaled in units of (pt/Qs) 2 , shows virtually no N y dependence 
in particular for the moderate pt region that dominates the integrated multiplicity. We 
may conclude that it has a weak dependence on N y (for fixed m/Q s ) addressing a concern 
raised in Ref. [jj7fl . The numerical single particle spectrum in the UV region is sensitive 
to the lattice spacing a, which, for a fixed value of Q s Ra, translates into a dependence 
on the size of the lattice Nt- As shown in fig. [3|, the spectrum approaches the continuum 
ultraviolet behavior of Qg/p T with increasing lattice sizes from Nt = 128 2 to Nt = 512 2 
lattices. 



N=512, g 2 n A =1.0 GeV, m/g^ fl =0.1 g 2 l i fl =1.0,m/g 2 n A =0.1 , N y =50 




' ' ' ' ' ' U ' ' ' ' ' 

2 4 6 8 10 5 10 15 20 25 

p T /Q s p T /Q s 



Figure 3: Left: The single particle spectrum (scaled by (pt/Q s ) 2 ) for different N y . Parameters 
are g 2 /j, — 1 GeV, Nt = 512, IR cutoff fixed at O.lg 2 ^. Right: The single particle spectrum (scaled 
by (pt/Qs) 4 ) for different Nt- The continuum spectrum has a Qt/Pr behavior at pt 3> Q s ', the 
convergence to this continuum behavior is clearly seen with increasing lattice size. 



We studied the m/Q s dependence of the single inclusive spectrum for a fixed large N y = 
50 and a 512 x 512 transverse lattice. The result is shown in fig. ^| (left) . For small m/Q s , the 
dependence of the spectrum on this quantity is quite weak-changing m/Q s by a factor of 2 
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shows virtually no change in the spectrum. However, some dependence is seen when m/Q s 
is increased, albeit the dependence is relatively weak in the pr ~ Q s region which gives the 
dominant contribution to the integrated multiplicity. It should be kept in mind that the 
dependence of the spectrum on m/Q s as shown is weaker than a logarithmic dependence 
even in the region where the dependence is the largest. It might seem counterintuitive that 
the value of an infrared scale m affects the gluon spectrum at such high momenta. One way 
to understand this is to keep in mind that the unitarity of the Wilson lines imposes a sum 
rule on their correlators, eqs. (fDJI) and (|3720|) . Namely, because C f (xt — Yt = 0) = N c , 
the integral of the momentum space correlator J ^?\htC f {k.x) = (^) 2 N C is a constant. 
Thus a modification of the distribution in the infrared will also affect the UV behavior, 
which is also reflected in the gluon spectrum. 

Finally, we plot in fig. Q (right) the dependence of the single particle spectrum on the 
other dimensionless combination of scales g 2 ^L for fixed m/g 2 fi. (This corresponds to a 
very good approximation to fixed m/Q s .). Virtually no dependence is seen on this quantity, 
confirming the expectation that the single inclusive spectrum is completely insensitive to 
physics on scales corresponding to the size of the system. In summary therefore, the single 
particle spectrum is most sensitive to physics at the scale Q s , weakly sensitive to physics 
on the scale m and completely insensitive to physics on the scale 1/R. 
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m/Q s =0.05 
m/Q s =0.1 
m/Q,=0.2 
m/Q s =0.3 
m/Q,=0.5 



© 



iiti 



4 6 

Pt'Qs 



10 



& 4 [■ 

2' 

« 3 
g 

s 2 



9 J y© 



CMa-1-0 ■ 

& A =2.0 • 
dl fl =3.0 

l>A=4.0 © 



3p_ 

©5i 



Pt'Qb 



Figure 4: Left: The single particle spectrum for different m for TVj, = 50, Nt = 512, and Q s = 
1 GeV. The spectrum is again scaled by (pr/Qs) 2 - Right: The dependence of the single particle 
spectrum on g 2 [i for a fixed raj g 2 /i. 



4.2 Double inclusive spectrum 

We now turn to the main focus of this study-the non-perturbative double inclusive gluon 
spectrum in A+A collisions. As in the single inclusive case studied in the previous section, 
we will first examine the sensitivity of the spectrum to the lattice parameters spelled out 
in section Because the sensitivity of the single inclusive spectrum to these quantities 
is known to be weak, we will directly plot the dependence of K2 (see eq. (|1.4| )) on some of 
these quantities. After exploring the sensitivity of our results to lattice artifacts, we will 
quote results for the physical value of K2, and discuss the implications of our results. 

In fig. H we plot the correlated part of the double inclusive spectrum as a function of 
Pt/Qs, the momentum of one of the tagged gluons, while holding the transverse momentum 
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N=256, Q s =1.0 GeV, N =50, 2.9<q T /Q s <3.1 



N=256, Q B =1 .0 GeV, N =50, 2.9<q T /Q s <3.1 



0.0014 
0.0012 
0.001 
0.0008 
0.0006 
0.0004 
0.0002 




♦ 

.1, 



m/Q s =0.05 i — B- 
m/Q s =0.1 : 
m/Q s =0.2 
m/Q s =0.3 i • 
m/Q s =0.4 
m/Q„=0.5 



« I i 




Pt'Qs iog(P T /Q s ) 

Figure 5: Sensitivity of the correlated part of the double inclusive spectrum to the value of IR 
cutoff m/Q s . Left: Linear plot scaled by (pt/Qs) to examine the anticipated perturbative behavior 
for large pr- Right: Log-log plot of double correlated part of the double inclusive spectrum. The 
results are plotted for a small bin in qx around 3 Q s . 
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Figure 6: Left: Dependence of K2 on the transverse lattice size Nt as a function of pr/Qs for 
fixed qr/Qs- Right: Dependence of «2 on the number of points used in the construction of the path 
ordered exponential, N y , as a function of pr/Q s for fixed qr/Qs- 



Qt/Qs of the other gluon fixed. The spectrum shown, plotted for different values of m/g 2 fi, 
is scaled by a factor (pt/Qs) 4 in fig. [| (left). From the perturbative computation discussed 
in section |!| we expect this scaled spectrum to go to a constant value at large pt/Qs- We 
note that it does so approximately, keeping in mind that the spectrum at large pt/Qs is 
especially sensitive to lattice artifacts. We also note that the spectrum is weakly dependent 
on m/g 2 fi. In fig. [5] (right), we also plot the correlated double inclusive gluon spectrum 
as a log- log plot. In addition to the power law tail, we observe a qualitative change to a 
softer px spectrum for px < 3Q S . 

4.3 Determining k% 

In fig. ||, we study the dependence of K2 on the number of transverse Nt x Nt lattice sites 
and on N y , the number of points in the longitudinal direction used to construct the path 
ordered Wilson lines. In the former case, hardly any dependence is seen thereby indicating 
a rapid convergence to the continuum limit. For longitudinal lattices, one observes a 
strong dependence for small N y , but a rapid convergence thereafter, with virtually no N y 
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dependence seen for N y > 10 onwards. The dependence of K2 on g 2 \iL for a fixed value of 
m/g 2 fj, is shown in fig. [7|. As anticipated, the dependence of K2 on g 2 [iL is rather weak. 
This further confirms, as suggested previously, that the physics in the infrared is, for a 
finite m, insensitive to the size of the system Ra- 

Now that we have established the convergence of K2 as a function of N y and Nx is 
robust, we should remind the reader that the continuum value of K2, being dimensionless, is 
most generally expressed as K2(px/Qs, qr/Qs, A<^, Q s Ra, m/Q s ). In general K2 is a function 
of two 2-dimensional vectors. When one takes into account rotational symmetry, it is a 
function of three variables and there are many ways to plot such a function. The general 
structure is illustrated in fig. ||, where we show K2 as a function of |pt — qr| and \px + q-r] • 
There is a delta-function peak at pt = qr which is clearly visible. The peak is left out in 
the right hand plot of fig. H to better illustrate the remaining structure. As mentioned in 



sec. 



PT 



3.1, there is also an enhancement in the correlation in the back-to-back configuration 
= — qr, which can be clearly seen in fig. & 
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Figure 7: k-2 for differing values of g 2 g.L for a 
fixed m/g 2 g, which is equivalent here to a fixed 
m/Qs- 



We now address the dependence of K2 
on m/Q s . We saw in fig. [| that the double 
inclusive spectrum had a weak dependence 
on m/g 2 fi, but the single inclusive spectrum 
had a slightly stronger dependence. The re- 
sulting effect on K2 is shown in fig. |9| as a 
function of px with qx averaged over the 
interval [2.9 Q s , 3.1 Q s ] and (IS iX function of 
the angle between px and qr- One observes 
that K2 decreases rapidly with increasing 
m/Q s but converges to K2 ~ 0.5 for larger 
m/Q s . The interpretation of the m/Q s de- 
pendence of the results will be discussed fur- 
ther in the next section. 

We can establish from the r.h.s. plot in fig. ^ that K2 is nearly constant as a function 
of Aip; the strength of the correlation is only weakly dependent on the relative azimuthal 
angle between pairs of gluons. This result confirms the conjecture in Ref. [22]. Turning 
to the dependence of K2 on px/Qs, Qt/Qs, the dependence is not entirely flat as assumed 
in Ref. [22]; nevertheless, despite some structure, the variation of K2 with pr/Qs is on the 
order of 20% for p T /Q s < 4. 

Figures 10 and |ll] show the effect of taking different configurations of the "trigger" 
and "associate" transverse momenta px and qx- In fig. [H] (left) the momentum qx is fixed 
to three different narrow bins and K2 is plotted as a function of px- Figure |l0| (right) shows 
the dependence of K2 on the angle between p^ and qr for three different combinations of 
Pt and qT', either both around Q s , both around 3Q S or one at Q s and the other one at 
3Q S . In fig. 11, we show K2 for px and qx having the same value, but with a finite angle 
between the vectors px and q-r; showing a strong increase in the correlation for smaller 
values of the momenta. 
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Figure 8: Three dimensional plot of k(|pt — qr|, |Pt + 1t\)- The plot on the left shows the delta 
function peak at = qr- On the plot on the right this peak has been left out by reducing the 
range (the |pt — q-rl/Qs-axis starts from 0.1) to show the structure around the away-side, i.e. close 
to pt = — qr- 



N T =256, Q s =1.0 GeV, N =50, 2.9<q T /Q s <3.1 
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Figure 9: k-i for different values of m/Q s with Q s = 1 GeV, at N y = 50. Left: K2 as a function 
of pt for qr in a small bin around 3<5 S - Right: K2 as a function of Atp with both pr and qr in a 
small bin around 3Q S . 



5. Discussion and physical interpre- 
tation of results 

We shall now discuss the physical implica- 
tions of our numerical results for the double 
inclusive gluon distribution in the Glasma. 



In the Glasma flux tube picture [22|, as 



noted in eq. (1.4) previously, it was conjec- 
tured that the correlated two gluon spec- 
trum C(p, q) can be expressed as 



C 2 (p,q) 



dN 
dy p d 2 p T 



dN 



dy q d 2 q T 



K2 



S±Qi 



(5-1) 
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Figure 11: k as function of momentum at equal 
momenta and specified angular separation. 



where k 2 is a non-perturbative constant which behaves parametrically as N c . The per- 
turbative computation of k 2 has a logarithmic dependence on pr/Qs and qr/Qs- Modulo 
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Figure 10: Left: k as a function of absolute value of one of the momenta, the other being fixed to 
the specified value. Right: Angular dependence of k at different absolute values of the momenta. 
Note absence of (5-function (point at A(p = outside plotted range on the y-axis) at unequal 
momenta. 



this logarithm, one estimates perturbatively 9 that K2 ~ 0.4. 

In comparing our non-perturbative results to the perturbative estimate, we first note 
that K2 is weakly dependent on Aip and pr/Qs, Qt/Qs, just as conjectured in Ref. p2| . 
Our result for K2 is larger than the perturbative estimate, but it has a sensitivity, as shown 
in fig. ||, to m/Q s . As noted previously, the double inclusive spectrum by itself has a 
relatively weak dependence on m/Q s and the dependence on m/Q s in the ratio K2 mostly 
comes from the single inclusive spectrum. Let us now discuss the physical interpretation 
of this dependence. 

Recall that in our computation m appeared as a cutoff inserted into eq. ( |3.9D to invert 
the Laplace operator. It regulates the long distance Coulomb tails of the color field 10 . 
In the context of the MV model, this scale appears as an additional external parameter, 
whose value is not determined within the model itself. Thus it is not completely clear 
whether the scale m should be thought of as the confinement scale or a scale of order Q s . 
In a straightforward implementation it is natural to think of m as a cutoff imposing color 



neutrality at the size of a nucleon. This is indeed the picture used for instance in Refs. [p~8 | . 

This picture does not, however, take into account color screening effects coming from 
quantum evolution effects (virtual and radiative corrections) that are described by the 
JIMWLK and BK equations. At RHIC energies they may not be very important but will 
be extremely important at LHC energies. The effect of color screening of the correlators 
has been considered previously Q] (see also the discussion in Ref. p3[) and it is argued 
that quantum evolution effects regulate the infrared behavior of the color charge density 
correlator at a scale that is also parametrically of the order of Q s . This is not a purely 
classical effect arising from nondinear Yang-Mills dynamics but from a combination of 
rescattering and quantum evolution effects. To determine the "correct" value of m in this 
context it is thus necessary to go beyond the MV model and include the effects of high 



9 On account of incorrect factors of 2 and it, the number quoted in Ref. |22| is an order of magnitude 
larger. The correct perturbative value (modulo logs) is quoted in Refs. ]3^, |32j], 
10 The propagator for the 2 dimensional Laplace operator is a logarithm. 
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energy evolution. A systematic numerical study of the two gluon correlation when high 
energy evolution effects are included is beyond the scope of this work. 

A practical way to address this question in our simulation is to note that for a geo- 
metrical flux tube picture to be valid the two particle cumulant should be given by the 
size of the flux tube Sf t (representing the transverse range of color correlations) divided by 
the transverse area of the system times a number of order unity. We can make this logic 
explicit by expressing the size of the typical flux tube as 

C 2 (p,q) 5ft 



dN \ / dN \ S± 
dy P d 2 p T / \ dy q d 2 q T , 



(5.2) 



where 5ft = k 2 (- • • , m /Qs)/Qs- For the range of m/Q s = 0.1-0.5 considered, k 2 ~ 2 - 0.5 
converging rapidly to the latter value for increasing m/Q s . The effective scale governing 
correlation strengths of unit strength is then 1/y/Sft ~ 0.7— 1.4 Q s . This scale, although 
not numerically very large, is nevertheless a semi-hard scale and corresponds to transverse 
distances much smaller than the nucleon size, thereby confirming the picture of early times 
in the Glasma as classical field configurations that are coherent over transverse distances 
much smaller than a nucleon. 

Let us then turn to a comparison with the RHIC data on two particle correlations. 
The experimentally observed quantity is Ap/ ^fp^(Aip); in our framework, for Aip = 0, it 
can be expressed as 

Ap . dN C 2 (p,g) ( 1, 
(Atp = 0) = — • — — r— — — r (7b ) 



Vp^- d y ( AI i> ) dN > ) ib' 

\ dy p d<*p T / \ dy q d^q T / 

= ^(7B-^), (5.3) 
where = k 2 /13.5 for an SU(3) gauge theory and 75 is the average radial boost in the 



framework of a blast wave model. From the RHIC data 29], one can estimate that 



Ap/y/p Te f_(Aip = 0) = 1/w27tcj2, with Gtp = 0.64. Combining this with eq. ( |5.3| ), one 



obtains 



for a s = 0.5 and where the superscript denotes that this is a crude estimate extracted from 
experiment using a blast-wave parametrization. For an average blast wave radial velocity 
V r = 0.6, this gives ftf w ~ 1.5; for V r = 0.7, one obtains ~ 1- There is considerable 
variation therefore in the value of w obtained from the ridge data due to the final state 
flow parametrization. 

An alternative method to extract k 2 is to compare the expression for the negative 



binomial multiplicity distribution (cf. eqs. (2.11 )— (2.12) in the Glasma [ 33 [ to PHENIX 



data [34] on the same. From this comparison, one obtains 

4 BD ~ 3.9. (5.5) 

One sees therefore considerable variation in the extraction of k 2 from experiment within 
the present framework. Within the many uncertainties, one can say at best it is a number 
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of order unity. Our study suggests that a coherent picture of such correlations from RHIC 
and higher LHC energies can in principle provide unique information on color screening of 
strong fields at early times in heavy ion collisions. 

An objective of this study is also to proceed in the opposite direction, namely, to 
determine K2 from a non-perturbative computation and use this as input into a dynamical 
space-time evolution model. As guide to this future program is the work in Ref. p9[ . 
where two particle correlations are extracted from the hydro dynamical evolution of flux 
tube structures in A+A collisions. 

6. Summary 

In this paper, we investigated the validity of the Glasma flux tube scenario of multi-particle 
correlations by performing a non-perturbative numerical computation of double inclusive 
gluon production in the Glasma. Our results were obtained by solving Yang- Mills equations 
on a space-time lattice for Gaussian distributed color source configurations with a variance 
proportional to the saturation scale Q^. Our results confirm key features of the Glasma 
flux tube scenario. Particles produced from coherent longitudinal electric and magnetic 
fields in transverse regions of size ~ 1.4/Q s -l/2 Q s , much smaller than the nucleon size, 
have correlations of unit strength. As in the asymptotic perturbative estimates, particles 
produced from the flux tubes are uncorrelated in the relative azimuthal angle between the 
particles, the observed azimuthal collimation being produced later from radial flow. The 
correlations show non-trivial structure in p? and qx, which smoothly go over to perturbative 
results in the limit of pt/Qs S> L Qt/Qs 3> 1- Our results for the two particle correlation 
strength are consistent with estimates of the strength extracted from model comparisons to 
RHIC data. A useful extension of this work is to incorporate our results into more detailed 
hydrodynamical models of the space-time evolution of initial state correlations. 
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A. Lattice formulation 

The Yang-Mills equations can be formulated as Hamilton's equations of motion. To pre- 
serve gauge invariance, they are solved numerically on a lattice where the degrees of freedom 
at a site i are the link variables 




(A.l) 
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where a is the lattice spacing on a transverse lattice. The appropriate discretized lattice 
(Kogut-Susskind) Hamiltonian in our case is given by 

aH = y / Tr E i E i + J-Re Tr U u )+ - Tr tt 2 + - V Tr U - , 

I r g 2 a \ N c J a r ; V 

(A.2) 

where Ei,Ui,ir and are dimensionless lattice fields that are matrices in color space, with 
E l = E a % t a etc. and the generators of the fundamental representation normalised as 
Tr (t a t b ) = 6^/2. 

The first two terms are the transverse electric and magnetic fields with the transverse 
plaquette defined as 

C/n(x T ) = C4(x T )f/ y (x T + e x )Ui{* T + e y )Ul^ T )- (A.3) 

The last two terms are the kinetic energy and covariant derivative of the rapidity component 
of the gauge field (j) = = —t 2 A v . The latter becomes an adjoint representation scalar 
following the assumption of boost invariance. For the parallel transported scalar field, we 
use the notation 

&(x r ) = 0i(x T )0(x T + ei)E//(x T ). (A.4) 

In the Hamiltonian, eq. (|A.2| ), there is a residual invariance under gauge transformations 
depending only on the transverse coordinates. 
The Hamiltonian equations of motion are 



(A.5) 
(A.6) 
(A.7) 



Ui = 


q 2 ■ 

i^-E l Ui 


(no sum over i) 


T 




<p = 


TIT, 




E x = 


— \u 


+ U x - y - h.c. 


E y = 


— \U 


+ Uy- X - h.c. 


7T = 




+ 4>-i - 2(j) . 



h.c. 1 — trace + 



T 

~ [ ~ 

T 



(Ai 



Gauss's law, conserved by the equations of motion, is given by 

J2 [uH*T - ei)^(x T - ei)Ui(xT ~ <*) - £?(x T )] - vr] = . 

i 

On the lattice the initial conditions (p.4|) are 



= Tr 



t a UUl + Uf)(l + Un- h.c 



E l = 0, 

= 0, 
-i 

±9^7 



(C/ i (x T )-l)([/f(x T )-C// 1 (x T ) 



+ ([//(x T - ei) - l) (t/?(x T - et) - Uh*T - ei)) - h.c. 



(A.9) 

(A.10) 
(A.ll) 

(A.12) 
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where U: ' in eq. (|A.9| ) are the link matrices corresponding to the color fields of the two 



nuclei {A\' 2 in eq. ( |3.5[ )). The link matrix Ui, which corresponding to the r > color field 
Ai, is determined by solving eq. 
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